#include <iostream>
#include <vector>
#include <fstream> // 包含 ifstream 所需的头文件
#include <complex> //使用复数的头文件
#include <chrono> //计算耗时
#include <../../include/fftw-3.3.10/fftw3.h>
#include "../../include/mfccs/taiic_mfcc_v3.h"


using namespace std;
using namespace Eigen;

#define M_PI 3.1416
const int winLength = fftSize;
const int num_filters = num_mfcc_coef * 2; 

//取对数与计算离散余弦变换（DCT）  
void dct(const MatrixXd& input,MatrixXd& output)
{
    int rws = input.rows();
    int cls = input.cols();

    MatrixXd tmp_out(rws, cls);

    double factor = M_PI / cls; 

    for(int m = 0; m < rws; ++m)
    {
        VectorXd temp = input.row(m);
        for (int k = 0; k < cls; ++k)
        {
            double sum = 0.0;
            for (int n = 0; n < cls; ++n)
            {
                sum += log(input(m,n) + 1e-10) * cos(factor * (n + 0.5) * k);//计算对数能量:log(input[n]+1e-6)
            }
            // 归一化  
            if (k == 0)
            {
                tmp_out(m,k) = sum * sqrt(1.0/cls);
            }
            else
            {
                tmp_out(m,k) = sum * sqrt(2.0/cls);
            }
        }
    }
    //取出前num_mfcc_coef特征数
    for (int i = 0; i < numFrames; i++)
    {
        for (int j = 0; j < num_mfcc_coef; j++)
        {
            output(i, j) = tmp_out(i,j);
        }
    }
}

MatrixXd melFilterBank(int fs, int fft_size, int num_filters)
{
    double lowFreq = 0.0;
    double highFreq = fs / 2;
    //计算梅尔滤波器组  
    double mel_min = (2595.0) * log10(1 + lowFreq / 700.0); // 0 Hz  
    double mel_max = (2595.0) * log10(1 + highFreq / 700.0); //Nyquist frequency  

    VectorXd mel_points(num_filters + 2);
    // 生成梅尔点  
    for (int i = 0; i < num_filters + 2; ++i)
    {
        mel_points(i) = mel_min + i * (mel_max - mel_min) / (num_filters + 1);
    }

    // 反梅尔频率转换  
    VectorXd hz_points(num_filters + 2);
    for (int i = 0; i <= num_filters + 1; ++i)
    {
        hz_points(i) = 700.0 * (pow(10, mel_points(i) / 2595.0) - 1);
    }

    VectorXi bin_points(num_filters + 2);
    for (int i = 0; i <= num_filters + 1; ++i)
    {
        bin_points(i) = static_cast<int>(floor((fft_size + 1) * hz_points(i) / fs));
    }
    //创建滤波器  
    MatrixXd mel_filter = MatrixXd::Zero(num_filters, fft_size / 2 + 1);
    for (int m = 0; m < num_filters; ++m)
    {
        for (int k = bin_points[m]; k < bin_points[m + 1]; ++k)
        {
            mel_filter(m, k) = static_cast<double>(k - bin_points[m]) / (bin_points[m + 1] - bin_points[m]);
        }
        for (int k = bin_points[m + 1]; k < bin_points[m + 2]; ++k)
        {
            mel_filter(m, k) = static_cast<double>(bin_points[m + 2] - k) / (bin_points[m + 2] - bin_points[m + 1]);
        }
    }
    return mel_filter;
}

MatrixXcd fftw(const MatrixXd& signal)
{
    int rws = signal.rows();
    int cls = signal.cols();

    MatrixXcd output(rws, cls);

    fftw_complex* in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * cls);
    fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * cls);
    for (int m = 0; m < rws; m++)
    {
        for (int n = 0; n < cls; n++) 
        {
            in[n][0] = signal(m,n);
            in[n][1] = 0.0;
        }

        fftw_plan plan = fftw_plan_dft_1d(cls, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
        fftw_execute(plan);

        for (int i = 0; i < cls; i++) 
        {
            output(m,i) = complex<double>(out[i][0], out[i][1]);
        }
    }
    return output;
}

//预加重
VectorXd pre_emphasis(const VectorXd& signal)
{
    double alpha = 0.97;
    VectorXd result(signal.size());
    result(0) = signal(0); //保持第一个样本不变
    for (size_t i = 1; i < signal.size(); ++i) 
    {
        result(i) = signal(i) - alpha * signal(i - 1);
    }
    return result;
}

VectorXd hanning_window(int length)
{
    VectorXd window(length);
    for (int n = 0; n < length; n++)
    {
        window(n) = 0.5 * (1 - cos(2 * M_PI * n / (length - 1)));
    }
    return window;
}

//归一化MFCC输出数据
void normalizeMFCC(const MatrixXd& mfcc, MatrixXd& mfcc_normalized)
{
    //计算最大值和最小值
    double minVal = mfcc.minCoeff();
    double maxVal = mfcc.maxCoeff();
    // printf("Max=%lf,Min=%lf\n", maxVal, minVal);

    //进行最大最小值归一化  
    MatrixXd norm_matrix = (mfcc.array() - minVal) / (maxVal - minVal + 1e-6);
    mfcc_normalized = (norm_matrix * 255.0).array().floor(); // 乘以255并向下取整  
}

// 计算 MFCC  
void compute_mfcc(const MatrixXd& signal,MatrixXd& out_uint8)
{
    //创建窗函数
    VectorXd window = hanning_window(winLength);
    //生成梅尔滤波器组 
    MatrixXd mel_filters = melFilterBank(sampleRate, fftSize, num_filters);
    //对信号进行加窗
    MatrixXd windowed_signal=signal.transpose();//480*400
    for (int i = 0; i < windowed_signal.cols(); i++)
    {
        windowed_signal.col(i) = windowed_signal.col(i).cwiseProduct(window);
    }

    MatrixXd fft_intput = windowed_signal.transpose();
    //fft变换
    MatrixXcd spectrum = fftw(fft_intput);   

    //计算能量谱
    MatrixXd power_spectrum(numFrames,fftSize / 2 + 1); //创建能量谱向量:400*241 
    for (size_t i = 0; i < numFrames; ++i)
    {
        for (size_t j = 0; j < (fftSize / 2 + 1); ++j)
        {
            power_spectrum(i, j) = abs(spectrum(i, j)) * abs(spectrum(i, j));
        }
    }

    //计算梅尔频率
    MatrixXd mel_spectrum(numFrames, num_filters);//创建梅尔频率向量矩阵
    for (int i = 0; i < numFrames; i++)
    {
        for (int j = 0; j < num_filters; j++)
        {
            mel_spectrum(i,j) = mel_filters.row(j) * power_spectrum.row(i).transpose();
        }  
    }
    //计算log对数+DCT变换
    MatrixXd out_dct(numFrames,num_mfcc_coef);
    dct(mel_spectrum,out_dct);

    //归一化到uint8
    normalizeMFCC(out_dct,out_uint8);
}


//读取音频数据->VectorXd
VectorXd loadAudioData(const char *filename)
{
    std::ifstream file_pcm(filename, std::ios::binary | std::ios::ate); //二进制形式打开，并且打开时文件指针立即置于文件末尾
    if (!file_pcm.is_open())
    {
        std::cerr << "无法打开文件 " << filename << std::endl;
        throw std::runtime_error("文件打开失败");
    }
    // 获取文件大小(单位:字节)，并转换为16位样本的数量
    streamsize size = file_pcm.tellg();
    if (size % 2 != 0)
    {
        throw std::runtime_error("文件大小不是16位样本的整数倍");
    }
    size_t nSamples = size / 2;
    file_pcm.seekg(0, std::ios::beg); //重置到文件开头

    //创建一个足够大的vector来存储数据
    std::vector<int16_t> pcmData(nSamples);

    //读取数据到vector
    file_pcm.read(reinterpret_cast<char*>(pcmData.data()), size);
    if (!file_pcm)
    {
        std::cerr << "读取文件时发生错误" << std::endl;
        throw std::runtime_error("文件读取错误");
    }
    //关闭文件
    file_pcm.close();

    VectorXd vec(nSamples);//vector<double> audio_f(nSamples,0.0);
    for (int i = 0; i < nSamples; i++)
    {
        vec(i) = static_cast<double>(pcmData[i]);
    }
    
    return vec;
}

void extract_mfcc(const VectorXd& audioD, MatrixXd& out)
{
    int numSamples = audioD.size();
    int TotalLength= (numFrames - 1) * frameShift + winLength;
    VectorXd vector_norm(TotalLength);
    vector_norm.setZero();//初始化为0

#if TAIIC_ASR  //自研算法进行归一化操作
    //计算均值与方差,零均值标准化,也称为Z-score标准化
    double mean_v = audioD.mean();

    double variance = (audioD.array() - mean_v).square().sum() / (audioD.size() - 1); // 使用无偏估计  
    double variance_v = std::sqrt(variance);

    //进行均值归一化和标准差归一化  
    VectorXd normalized_audioD = (audioD.array() - mean_v) / (variance_v + 1e-6);
    vector_norm.head(numSamples) = normalized_audioD;
 
#elif ACOUSTICS_ASR //声学所算法不进行归一化操作
    vector_norm.head(numSamples) = audioD;
#endif

    VectorXd pre_emphasized_signal = pre_emphasis(vector_norm);
    //1.分帧
    MatrixXd frames(numFrames,winLength);
    for (int frame_idx = 0; frame_idx < numFrames; frame_idx++)
    {
        /*printf("frame_idx=%d\n",frame_idx);*/
        int start = frame_idx * frameShift;
        int end = start + winLength;
        VectorXd frameData = pre_emphasized_signal.segment(start,end-start);
        frames.row(frame_idx) = frameData;
    }//for frame_idx

    //2.计算MFCC
    compute_mfcc(frames, out);
}







